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transformations 
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This q)plication claims priority benefit to U.S. provisional application number 
5 60/262,319 filed January 17, 2001. 

FIELD OF THE INVENTION 

The present invention generally relates to systems and methods for performing 
simulations. More specifically, the invention relates to systems and methods that 
perform calculations in multiple domains and use domain transformations to move 
10 between the domains. The transformations preserves a basic, underlying principle 
common to all domains, which is usefiil in simulating one or more properties of a 
multi-component fluid contained in a physical system. 

BACKGROUND OF THE INVENTION 

15 

Numerical simulation is widely used in industrial fields as a method of 
simulating a physical system by using a computer In most cases, there is a desire to 
model the transport processes occurring in the physical system. What is being 
transported is typically mass, energy, momentum, or some combination thereof. By 

20 using numerical simulations, it is possible to model and observe a physical 
phenomenon or several coupled phenomena. 

Basin simulation is of great interest because it provides quantitative prediction 
of the timing, volumetrics, quality and distribution of oil and gas in a petroleum 
system as well as pressure. From an exploration point of view, it provides a means to 

25 quantify the risk of drilling opportunities by prediction of the most probable volumes 
of petroleum accumulated, quality of petroleum, timing (petroleum generation relative 
to trap formation), targets of commercial petroleum accumulation based on sensitivity 
analysis of poorly constrained basin parameters and pressure. All simulations are 
constrained by interpreted seismic data, chrono-stratigraphy, geochemical and well 

30 log data. 
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An integrated basin simulator is complex. It is able to perform numerical 
modeling through geologic times of coupled basin processes such as sedimentary 
deposition, burial and compaction, structural evolution, uplift and erosion, thermal 
evolution, petroleum generation and expulsion, fluid migration, petroleum 

S accumulation and pressure estimation at the present day. One key feature is the 
presence of organic chemical reactions modeling the solid karogen decomposition 
coupled with temperature evolution that generates solid and flowing petroleum 
products. These petroleum products may also undergo fiirther cracking reactions as a 
fimction of temperature during fluid migration. Accordingly, the organic chemical 

10 reactions are coupled with mass transfer and fluid flow process. Large scale rock 
deformation associated with burial is also coupled with fluid flow processes. The 
pore fluid pressure may have a large effect on rock deformation and porosity and, at 
any instant, the molar vohimes of the flowing fluid phases stored in a rock are 
strongly influenced by porosity. Because these coupled processes are so complex, 

15 basm simulations can only be done using computers. The objective of basin 

simulation is to understand the complex interaction of these processes occurring in a 
basin to help minimize the risk of petroleum exploration and drilling. 

Modeling of these complex processes involves use of partial differential 
equations. As a means for numerically solving such equations, there are known the 

20 finite element method, the finite difference method, the finite volume method and the 
like. Regardless of which method is used, the physical system (i.e., the basin or the 
reservoir) to be modeled is divided into cells (a set of which is called a grid or mesh), 
and the state variables that vary in space throughout the model are represented by sets 
of values for each cell. In order to analyze the evolution of the basin that changes in 

25 time, it is necessary to calculate physical quantities at discrete intervals of time called 
timesteps, irrespective of the continuously changing conditions as a fimction of time. 
Basin simulation proceeds m a sequence of timesteps. Althougfh coupled, some of the 
processes can be evaluated sequratially within a global timestep. This results in 
lagged variables for some processes. For example, thermal calculation may be 

30 performed prior to the resolution of rock motion coupled with fluid flow that can then 
be performed under isothermal conditions. 



wo 02/057901 



3 



PCT/US02/01041 



Simulations based on immiscible fluid flow, for which there is no component 
exchange among the fluid phases (aqueous, liquid petroleum (oil) and vapor 
petroleum (gas)) are often not as representative of multi-component, multi-phase fluid 
flow as compositional simulations. Compositional simulations based on 

5 compositional fluid flow are well known to reservoir simulation as a process of 
inferring the behavior of a real reservoir from the performance of a model of that 
reservoir. The compositional model describes reservoir hydrocarbon as a multiple- 
component mixture. Gas and oil phase properties and equilibrium are calculated from 
pressure and composition dependent correlations or, more typically, from suitable 

10 equations of state (EOS). Several EOSs have been developed and are in use today, 
including for example the Redlich-Kwong EOS and the Peng-Robinson EOS. 

Compositional reservoir simulators using an EOS to describe the phase 
behavior of multi-component fluid mixtures are expensive to use because of the large 
number of iterative phase equilibrium calculations and large computer storage space 

15 required. The number of equations having to be solved in EOS calculations is 
proportional to the number of components in the fluid. To limit the computational 
time of compositional reservoir simulation, a common practice is to pseudoize the 
fluid description. In the pseudoization, the pure compounds are grouped into a 
number of component groups, termed pswdocomponents. The pseudocomponents 

20 are treated as if they were pure components in subsequent reservoir simulations. 

Basin processes are relatively slow so that the assumptions and constraints of 
phase equilibrium based on EOS equations^ as mooned above, equally apply to 
basin simulatioa Afiscible compositional fluid flow is important to basin simulation. 
The exchange of the components between the phases, which depends on pressure, 

25 temperature and composition, can afifect numerous phenomena that cannot be 

observed with an immiscible flow model. When light (gas) components are present in 
the liquid petroleum phase, the viscosity and the density of the phase are decreased. 
Consequently, the liquid petroleum phase migrates faster than in the immiscible 
situatiort On the other hand, the presence of oil components in the vapor petroleum 

30 phase will tend to slow down the migration of the phase. Changes in pressure 

conditions during migration can lead to the exsolution of the light components that 
migrated within the liquid petroleum phase. Overall, the timing of migration, the 
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number of phases during migration, and the distribution and composition of the 
accumulations may be affected. Lnmiscibility of the aqueous phase may also not be a 
reasonable assumptioa Light petroleum components, in particular methane, are 
soluble in water or brine with decreasing solubility as salinity increases. Modeling, 
5 detecting and predicting the presence of low saturation gas and/or 

'fizz' (gaseous) water are important e?cploration problems in determming commercial 
gas accumulations. 

In the porous medium^ component mass balance and momentum balance based 
on generalized Darcy's law are used for fluid flow process. Persons skilled in the art 

10 of multiphase fluid flow recognize that the generalized Darcy's law requires 
knowledge of the number and properties of the flowing phases, such as intrinsic 
density, viscosity, capillary pressure, mole fraction and composition in terms of the 
flowing components. The compositional property is a key parameter, which is 
typically expressed in the form of a non-square matrix. With immiscible fluid flow, 

15 phases and components are the same, and their properties can be evaluated using 

equations based on observations. With compositional miscible fluid flow, however, it 
is desirable to take into accoimt the exchange of components among the phases based 
on phase equilibrium and flash calculations. Phase equilibrium and flash calculations 
can also provide information about the physical properties of the phases which are 

20 required to perform fluid flow calculations. 

Several kerogen decomposition models have been proposed that use generic- 
generated flowing petroleum products, such as gas and oil, or dry gas, wet gas and oil, 
or dry gas^ wet gas, light oil and heavy oil. These models depend on the region of 
interest, the type of kerogen present, the availability of calibrated data, and the 

25 method of characterization. The flowing components can undergo secondary 
cracking as a function of temperature. Products generated by the kerogen 
decomposition and their secondary cracking provide the source terms for the fluid 
flow. These flowing components can be defined in terms of EOS components or 
pseudocomponents by mappings that assure component mass balance. The mapping 

30 is based on detailed analysis of petroleum products known to be generated by similar 
kerogens and must be compatible with some basic parameters used in the kerogen 
model. The number of EOS components is typically greater than or equal to the 
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number of flowing kinetic components. For example, the wet gas can be defined as a 
composition (based on mole fraction) of methane, ethane, propane, butane and 
pentane. With this set of conditions, compositional fluid flow process will involve 
transformations from one domain of components to a second domain of components, 

5 and reciprocally. Such direct and inverse transformations have to be peifi)rmed in 
such a way that mass balance is preserved overall and that the non-square 
composition matrix of the phases in terms of the flowing components can be defined 
with its inherent properties. The increased computational cost of computing 
transformations to move between the domains may at first seem prohibitive, but the 

10 transformational cost may be more than ofi&et by computational savings that result 
from reduced domain size. 

A need exists for a method of performing miscible conqsositional fluid flow 
with a set of flowing components different from, but related to, a set of EOS 
components used for evaluating the exchange of the components among the phases. 

15 The two sets of components are different in terms of number and are assumed to be 
linked by a mapping, constant in time and space, that preserves mass. The phases are 
assumed to be identical in both domains. The heart of the problem is to evaluate the 
composition of the phases obtained in the EOS domain in terms of the flowing 
components with all its inherent properties. 

20 It is noted, however, that the transformations linking the various domains are a 

potential source of numerical error and instability. It is therefore desirable to provide 
a robust transformation method that enables one to build a practical^ coherent 
simulation of a complex system. 

25 SUMMARY OF THE INVENTION 

This invention relates to a method for evaluating the composition of the 
flowing phases in terms of a set of flowing components that are different from, but 
related to, a set of EOS components or pseudocomponents used to obtsdn the number 
and properties of the flowing phases as a fiinction of pressure and temperature based 
30 on phase equilibrium and flash calculations. The flowing components, which can be 
reactive in their own domain, are related to the EOS components and/or 
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pseudocomponents through a mapping, based for example on mole fraction, that 
preserves mass and is constant in time and space. The number of flowing components 
is less than or equal to the number of EOS components. One or several flowing 
components can be an EOS component and an EOS component can be present in 

5 more than one flowing component. The method is based on the definition of a non- 
square matrix that expresses the distribution of the components among the phases and 
can be evaluated ia both domains. Other matrices can also be defined that can be 
used for other applications. 

While much of the ensuing description is directed toward oil field basm 

10 simulation, it is understood that the disclosed method will also apply to reservoir 
simulation. 

DESCRIPTION OF THE DRAWINGS 

A better understanding of the present invention can be obtained when the 
following detailed description of the preferred embodiment is considered in 
conjuction with the following drawings, in which: 

Fig. 1 is a block diagram of a preferred fluid flow simulator embodiment. 
Fig. 2 is a diagram of a computer that may be configured to effectuate 
simulations according to the present invention. 

DETAILED DESCRIPTION 

An improved method is herein described for estimating one or more properties 
25 of a multi-component fluid contained in a physical system containing two or more 
phases. In the first step, the physical system is equated in at least one dimension to a 
muhiplicity of cells. For each cell, the multi-component fluid is characterized using a 
property of a first set of con^onents and the characterization is represented by a first 
vector. For each cell, the first vector is transformed to a second vector using a fkst 
30 matrix (first transformation matrix), the first matrix being indicative of the 

distribution of the first set of components and the second vector being representative 
of a property of a second set of conq)onents greater in number than the first set of 
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components. The second vectors are then used to determine the number and 
properties of the phases present in each cell. Elements of a second matrix are then 
determined that expresses distribution of the second vectors among the phases. The 
elements of a third matrix are then determined that expresses distribution of the first 

S vectors among the phases. A fourth matrix is then determined that expresses the 
composition of the phases in terms of the first vectors. The fi}urth matrix is then used 
to perform fluid flow calculations to estimate one or more properties of the multi- 
conqponent fluid. In practicing the invention, certain operations are pefi)nned in one 
domain and other operations are performed in a second domain. The first domain has 

10 a set of components that characterize the muhi-component fluid and the second 
domain has a second set of components, greater in number than the first set, that 
characterize the multi-component fluid. It is assumed that the phases in both domains 
are the same. 

The following description makes use of several mathematical symbols, many 
15 of which are defined as they occur in this description. Additionally, for purpose of 
completeness, a Ust of symbols containing definitions of symbols used herein is 
presented in an Appendix which is located at the end of this detailed description. In 
order to differentiate the variables defined in both domains, the superscript ^'^^is used 

to indicate the flowing components domain while the superscript ^'^is used for the 
20 EOS components domain. When no superscript is present, it is implied that the 
relations can be defined in either domain. An indicial notation is maintained to 
differentiate the variables related to components or phases. Lower-case, Greek-letter 
subscripts indicate a component while upper-case, Latin-letter subscripts indicate a 
phase. Einstein's convention (i.e., sum over repeated indices) is not used. 
25 Summation over indices will be indicated by the symbol S with the lower and upper 
bounds given at the bottom and top of the symbol, respectively. Vectors are indicated 
by underlined letters, while matrices are indicated by double underlined letters. The 
usual matrix product is indicated by the symbol " • 
Definitions 

30 Before proceeding with the detailed description of the invention, a definition 

of the phase mole fiactions of components of the multi-component fluid (or 
distribution of the components among the phases) is first presaited. Consider a fluid 
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mixture made of fluid components distributed over Np fluid phases that 
corresponds to a total mixture molar density m , which is expressed as the total 
number of moles per unit basm/reservoir volume. It is understood that the number of 
components is greater than or equal to the number of phases, that is^ N^>Np. Let 
5 be the mixture molar density of component a and rjj be the mixture molar 
density of phase J . By definition of the mass balance, the total mixture molar 
density is equal to the sum over all the molar densities of the components and to the 
sum over all the molar densities of the phases: 

i^c AT, 
«=l J=l 

10 It is noted that equation (1) can be written with any unit. The mixture molar density 
(m„) of component a can be defined as the sum over the phases of the mixture 
molar densities of component a in phase/, denoted , that corresponds to the total 
number of moles of component a in phase y per basin/reservoir volume. Similarly, 
the mixture molar density (tjj ) of phase/ can be defined as the sum over the 

15 components of the mixture molar densities of component a in phase J (d^). The 
relations are then: 

'«a=i;^«/>Va€[l,JVj; and nj=^Zd^>'^J 4^Np] (2) 

The mole firaction of component a in phase J , noted , is the ratio of the mixture 
molar density of components a in phase J {d^ ) to the mixture molar density of 
20 phase J^ftj). As such, is a numb^ equal to zero or positive but less than or 
equal to 1 . For any phase J , the sum over the components of the mole fi:actions 
is 1 as derived using equation (2): 

^a^=— ,0^x„,^l; with j;,x^ = ^ =^^^lVJe[lM (3) 

The mixture molar density of components a, m^^ can then be expressed as a function 
25 of the mixture molar density of phase J, Tfj , using equation (2) as: 
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Such a relation can be expressed in a matrix form. Let ^ be the iV;, vector \^ose 
5 elements are the jw„ , 7 be the vector'wAose elements are the t/j , and x be the 
non-square JV^ xNp matrix whose elements are the . Equation (4) can then be 
written as: 

m=xi7. (5) 
The phase (J) mole fraction of the component a, noted yj^, can be defined as the 
10 ratio of the mixture molar density of component a in phase 7 (J^ ) to the mixture 
molar density of component a {m„ ) . As such, is a number equal to zero or 
positive but less than or equal to 1. For any component a , the sum over the phases of 
the mole fractions yj^ is equal to 1 as derived using equation (2): 

15 The mixture molar phase densities % can then be related to the mixture molar 

densities of the components by a relation dual of that e}q)ressed by equation (4) 
as: 



a 



which can be written in a matrix form similar to equation (S) as: 

\^ere y is the Np xJV^ matrix whose elements are the mole fractions yj^ . 

A relation betwem the mole fractions and yj^ can be defined using their 

definition given by equations (3) and (6), respectively, as, for any component a and 
any phase J : 



wo 02/057901 



10 



PCTAJS02/01041 




VaE[;iVj and \/J^[hNp] 



Va€[l,JV,] and Vy6[l,^p] 



(10) 



(9) 



Such relations allow, with the knowledge of the elements of the m vector and the tj 
vector, to compute the elements of the y matrix from those of the x matrix and 

5 reciprocally. These definitions apply to any field where the concept of moles can be 
defined. Persons skilled in the art would recognize that the elements of the vectors m 
and 1/ can be defined in many different ways as long as their elements verify a 
relation similar to equation (1). This would correspond to normalization of the 
elements of each vector by a given number These relations are independent firom the 
10 calculations that provide the x matrix or the y matrix. 

Persons skilled in the art of phase equilibrium and flash and petroleum 
reservoir science are fiuniliar with equations (1) to (5). Vector m is usually referred 
to as the feed whose elements are defined as the global mole fiactions of the 
components, that is, the total number of moles of component a divided by the total 
15 number of moles in the mixture. The vector tj is referred to as the phase mole 

fi-action, its elements bemg defined as the global mole firactions of the phase, Le., the 
total number of moles of phase 7 divided by the total numb^ of moles in the mixture. 
For a given pressure^ temperature and a set of components with their properties 
related to a given equation of state (EOS), phase equilibrium computations, based on 

20 the equalities of the fiigacties of the components in the phases, provide the number of 
phases and their properties among which are the mixture phase molar density 
(equivalent to the phase volume firaction) ijj and the elements x^ of the non-square 

X Np matrix (see the recent book by A. Firoozabadi, Thermodynamics of 
Hydrocarbon Reservoirs'", McGraw-Hill, 1999). Also provided are derivatives with 

25 respect to the mbcture molar density of the EOS components and pressure. Such 
derivatives are required to p^orm multiphase compositional fluid flow based on 
Darcy's law. The calculations are typically performed with only petroleum 
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components. Interaction of the water component with the petroleum phase and of the 
petroleum components with the aqueous/brine phase can also be considered (see the 
article of Li Y-K and Nghiem L. X, 1986, ''Phase Equilibria of Oil, Gas and 
Water/brine Mixtures from a Cubic Equation of State and Henry *s Law**, The 
5 Canadian Journal of Chemical Engineering, Vol. 64, p. 486-495). 

The mventors have discov^ed that the non-square y matrix or its elements 

yj^ can be useful in many different applications. The non-square Np x N„ matrix y 

whose elements are defined by equation (6), is not the generalized inverse of the non- 
square matrix x whose elements are defined by equation (3). A form of the pseudo- 

10 inverse of the matrix x is: 



where the superscript ^ mdicates the matrix transpose and the superscript indicates 
the inverse. With such a definition, assuming that the Np x Np square matrix 

(x^-x) is invertible, the NpxNp matrixproduct C* = y*x is equal to the identity 




15 matrix, and the non-square matrix Np x N^ verifies a relation of the form of 

equation (8), that is, y • m = 77 . However, there is no guarantee that the elements yj^ 

verify the properties outlined in equation (6) as required by the physics of the 
particular appUcation of the invention and of the very definition of the matrix. 
Therefore, although mathematically correct, for most applications use of the matrix 
20 y* for is not desirable. 

The non-square matrix y defined by equations (6) and (8) is mathematically 

sound and can be used to build useful matrices. Consider, for example, theiV^ x^^ 
square matrix C defined as the product of the>^ matrix times the x matrix, i.e., 




(11) 




(12) 
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The matrix C is not equal to the x Np identity matrix but, by construction, ^ is 

one of its eigenvectors, associated with the eigenvalue 1. Indeed, equations (S) and 
(8) yield: 

^ = 2^.^ = 7- x- ^ o 2 = £'2- (13) 

5 Although not equal to the identity matrix, the NpXNp C mattrix exhibit useftd 

properties. Using equations (12), (6) and (3) the sum of the elements of any of the 
column vectors of the C matrix is proven to be equal to 1, considering any phase K, : 

Up > i^j^ f \ f Up ^ u^ 

J= I J=l^o=l J af5 IV^JB I J OSS I 

Because of the condition that the number of phases is always equal to or less than the 
10 number of components (i.e., Np ^ N^X the square matrix C is non-singular and can 

thus be inverted. It can then be shown that the sum of the elements of any of the 
cohimn vector of the JV^ x Np C"^ matrix, inverse of the C matrix, is equal to 1. Let 



be the Kronecker symbol, that is, = 
theiSr^ X Np identity matrix. By definition of an inverse matrix: 

15 t <=> Z^uC;^- 



0, fori^ K * 



Equation (14) then yields for any phase K : 

Up Up Up ( ^? \ i^p 

ZZQ/<^=Z<yir=i = Z S^^^ c-A = Z(^A (15) 

1= 1J» 1 Lai J= l^ia 1 y 7« 1 

which is the desired result. Persons skilled in the art of numerical simulation will find 
applications for the C and C'^ square N^, x Np matrices, especially m problems 
20 related to mversioa 

TheiV^ X N^ square matrix£ is defined as the product of the x matrix times 

thej; matrix, that is: 

£= x y, o c,^= Zi'^ajyjfi' (16) 

^ 1 
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Such a matrix is of lesser interest than its dual counterparts C because for most 

basin/reservoir simulations the number of components will be greater than the number 
of phases (that is, N„ >Np). Such a condition implies that the matrix c will be 

singular or near singular (i.e., the ratio of the smallest eigenvalue to the largest 
5 eigenvalue will be less than the precision required for the calculations). Therefore, 
the matrix cannot be inverted and the properties relative to its inverse (as described by 
equation (15) for the inverse of the C matrix) do not apply in the components 

domain. Howev^, by construction, the vector^ is an eigenvector of the c matrix 
associated with the eigenvalue one (as seen from equations (5) and (8)): 
10 ?! = 2'2'^i*2!'— <=> ^ = £-^- (17) 

Second, the sum of the elements of any of the colimin vectors of thee matrix is equal 
tol: 

p. 



15 Compositional miscible fluid flow 

In simulations of miscible compositional fluid flow using generalized Darcy's 
law, the principal equations to be solved involve the balance of component mass and 
the balance of momentum. In the following, the equations are written in terms of the 
flowing components. For a given flowing component a^^^ , the mass balance is: 

f^^div[»iai^"")]=^ (19) 

where / is the time and div indicates the divergence operator taken with respect to the 

Eulerian coordinates, represents the vector of the particle velocity of the fluid 

component a^^^ and q^jl represents the source molar density rate relative to the same 

component. For basin simulation, the source t^m comes from the Idnetic model 
25 (decomposition of the kerogen in the source cells and possible secondary cracking in 
every cell), while for reservoir simulation, the source (or sink) term comes from the 
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preseace of injection or producing wells. The particle fluid velocity of the fluid 
compoment a^^hs related to the particle velocity V^^ of the fluid phase J as follows: 

•^^'"^'t'^Mt''- (20) 
The momentum balance for the fluid phase J results from Darcy's law as: 
5 FW.K(^)_^.g^ iP-P.'P,gz) (21) 

where V^^^ is the solid particle velocity which is zero in reservoir simulation; grad 
indicates the gradient operation with respect to the Eulerian coordinates; K is the 
absolute permeability tensor; (p is the porosity of the material; k^j is the relative 
permeability; Sj is the saturation (that is, the volume fraction); fij is the viscosity of 

10 fluid phase J ; P is the reference pressure (generally that of the aqueous fluid phase) 
while is the capillary pressure; pj is the intrinsic mass density of the fluid phase 
J ; ^ is the acceleration of gravity; and z is the depth. Introducing equation (2 1) into 
equation (19) using equation (20) and the relation between the mixture phase molar 
density of the phase J (tJj ) and the intrinsic phase molar density of phase / (4> ) , 

IS that is: 

nj^pSj4j (22) 
leads to the component balance equations to be solved for each mbcture molar density 
m^fj^ of the flowing component a^^^ : 

20 The resulting system of equations is highly non-linear and requires knowledge 

of the number and the properties of the phases as a function of the mixture molar 
densities m^J)^ of the chosen flowing components a^^\ the pressure P and the 

temperature T which is assumed constant during the time step. The resolution also 
requires additional equations that act as constramts relative to the saturation in 
25 reservoir simulation and to the porosity in basm simulation. Models for the variation 



dt 



-div 
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of the relative permeabilities and capillary pressures also need to be defined. The 
resolution can be performed either explicitly (IMPES type) or fully implicit (see, for 
example, the book by Aziz K. and Settari, A., 1979, Petroleum Reservoir Simulation^ 
Applied Science Publishers). Either way, it involves an iterative method that requires 
5 the evaluation of the so-called Jacobian matrix based on the derivatives of the phase 
properties with respect to the principal imknowns which m basin simulation are the 
mixture molar densities of the flowing components, the pressure P , and the 

porosity ip . Equation (23) illustrates that once a domain is chosen for the flowing 
components, the system can be solved only if the composition of the phases with 

10 respect to these flowing components represented by the elements jk^{^^ of the non- 
square matrix x^^^ are known, as well as the phase properties. 

In conventional reservoir simulation, the flowing components are chosen to be 
in the EOS domain which is also the domain in which the number and the 
properties of the flowing phases are evaluated. The present invention allows use of 

15 another domain in which the components are related to the EOS components by a 
transformation, constant in space and time, that preserves mass. In basin simulation, a 
natural definition of the flowing components is the flowing components generated by 
the decomposition of the kerogen in the source rock, which are also subject to 
secondary cracking as a function of temperature. A different domain could be chosen 

20 for when performing reservoir simulation. Using two domains in peforming a 
simulation decreases the size of the system to be solved (the number of flowing 
components being less than the numb^ of EOS components). 

Fig. 1 schematically illustrates the principal steps in carrying out an oil basin 
simulation in accordance with the practice of the present invention. A general 

25 overview is first described, followed by a more detailed analysis. 

One cycle of the loop depicted in Fig. 1 r^resents a cycle of a single time 
step. The choice of the kerogen decomposition model including secondary cracking 
defines the source terms (represented by block 90) within each cell of the domam at 
each time step i^ch play the role of initial conditions to the fluid flow simulation at 
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each time step. The number of flowing components in each cell is then determined by 
solving the fluid flow equations. 

As indicated in block 100, the definition of the source terms 90 implies a 
choice of the nature and number of the flowing components as represented by the 

5 mixture molar densities casted in the vector m^^^ It is also implicitly assumes that 
these flowing components have been defined in terms of a certain number of EOS 
components. In block 102, the flowing components are transformed mto a set of 

predefined EOS components as represented by the mixture molar densities vector nl''^ 
in block 103. This transformation is by carried out by a transformation matrix 7^^^ 
10 of the type described by equation (24) below. In block 1 04, an EOS, such as Peng- 
"Robinson model, is applied to each EOS component to determme the phase 
distribution equilibrium of each EOS component in each of the cells. This determines 

the number of phases, Np ; the mixture phase molar densities rf*^\ and the EOS phase 

composition matrix x^'^ . Block 106 represents the resulting EOS component phases 
15 that results fi-om the calculations of block 104. These phases are preferably expressed 
in the form of EOS component fi-action matrix (equation (10)). 

To reduce simulation complexity, the remaining simulation operations are 
performed on flowing components. To restate the EOS component fi:action y m terms 

of the flowing components, an inverse transform (block 108) is applied to the 

20 matnx of block 106 to determine the phases in terms of flowing components (block 
1 10) as represented by the flowing component firaction matrix of block 108. 

Assuming that the phases are the same in both domains, the inverse transform takes 
the form of equations (32) as given below. 

In block 112, operations are performed on the flowing component phases of 
25 block 110 to determine a new distribution of flowing components. The operations in 
block 1 12 preferably include determining partial derivatives of key variables with 
respect to unknowns and a solution to fluid flow equations, using contributions firom 
source rocks and molecular cracking phase extraction, and other suitable operations. 
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Referring now to Fig. 2, the simulation steps outlined in Fig. 1 are preferably 
implemented as a software program executed by a computer 202. Simxilation results 
may be displayed on a monitor 204, and interaction with a user is accomplished via 
monitor 204 and an us^ mput device 206. The software may be stored on computer 

5 readable information storage media 208, or accessed over a network or other suitable\ 
in&rmation transmission medium (not shown). One of skill in the art will recognize 
that computer 202 incliides a processor and a memory, and that the processor may 
retrieve the software from memory. Execution of the software configures the 
processor to opiate on digital signals stored in memory and manipulate those signals 

10 so as to carry out the simulation in the manner specified by the software. 
Detailed analysis 

The principal steps in carrying out the present invention as illustrated in the 
flow chart of Fig. 1 are now presented in detail. 

At block 100, it is assumed that a choice has been made for a set of 

15 JV^^^flowing components different from the number of EOS components. Implicitly, 

the pressure and temperature are considered as known and the choice of equation of 
state has also being made. The flowing components have also been defined in terms 

of JV^'^ EOS components, for example, based on mole firactions. The definition is 
assumed to be constant in space and time. It is understood that the number N^^^ of 

20 EOS components is greater than or equal to the number flowing components. 

It is also imderstood that an EOS component can be present in several of the flowing 
components and/or that a component can be both an EOS component and a flowing 
component. The transformation from the flowing components domain to the EOS 

domain is defined through a non-square -A^^'^x N^^ matrix (7^"^^ with elements of the 
25 form T^^j/) ) that relates the mixture EOS component molar densities th^)) , casted in 
the J\^') vector n^'\ to the mixture molar densities nJj}^ of the flowing components^ 
casted in the iV^^^ vector n^^\ as: 
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Such an operation is perfonned in block 102 with the result casted in block 103'. 

The elements T^^l^o are mole fractions so that they vary between 0 and 1, 
bounds included, and expressed as: 

O^Tjfi^^l. (25) 

S Theronservationof massmbothdomamcanbee}q)ressedasthee^^ 
over the component molar densities in both domams, e^)ressed as. : 

Z "^1= Z "iSi (26) 

that requires for the sum over the elements of any of the columns of the 7^^^^ matrix to 
10 be equal to 1: 

2 ^Sc/) = 1- (27) 



The transformation matrix r^^^ is assumed constant in time and space. 

At block 104, phase equilibrium calculations are performed with the EOS 
components for the given pressure and temperature known in block 100. These 
IS calculations provide the number (iVp) of the phases as well as some of their 

properties (intrinsic and mixture) among which are the mixture phase molar densities 
Tfji^ that can be represented in a iV^ vector rf"*^ and the composition of the phases 

represented by the matrix x^^^ whose elements x^^^ verify the relations given in 

equation (3). The component and phase mixture molar density vectors rj^'^ and 
20 are related through equation (5). 

At block 106, in the EOS domain, the elements y^J^^y of the non-square 

Np X JV^*^ matrix are calculated using equation (10), i.e., y^J^^^^ 'fj^ / • 
By construction, these elements verify the properties given in equation (6), i.e, 

0^ y% <. 1, VaWe[l,iV<;)] and Vye[l,JV^] (28) 

25 and Z&=l>^-^^^^p]- (29) 
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The mixture phase and EOS component molar densities are then related by an 
equation of the type of equation (8), expressed as: 

2^-)=/).mW. (30) 

The assumption of equality of the phases applies to the mixture phase molar densities 
5 rf''^ and rf^"^ . Using the transformation between the component molar densities 
(equation (24)) and the relation (30) written in both domains, a relation is fo\md 
between the non-square Np x matrix and the non-square Np x 1^/^ matrix 

10 By identification of the coefiGicients of the elements m^on both sides of the above 
equation, a particular solution is obtained as: 

= o;;J2,= 2 ^i^o^ico (32) 

Sudi an operation is performed in block 108. 

Because all the elements T^)jj) and j^^^j^j are numbers between 0 and 1, 

15 bounds included, as described by equations (25) and (28), the elements j^f^^ are also 

numbers between 0 and 1, bounds included. Also, the sum of all the elements of any 
column of the y^^^ matrix is equal to 1 as can be proven using equation (27): 

i*.= I fiy^,l3a»= f 2i^=i.vye[i.^,] (33) 

The matrix >;(^Mefined by equation (32) has then all the required properties and can 
20 be used to relate the phase and component molar densities in the kerogen domain as: 

2(/)=y/).y/) (34) 

In block 1 10, the elements j^j)^^ of the non-square N^^^ x Np composition 
matrix x^-^^ are calculated using equation (9), that is, jc^^ = y^J^ jtij. By 
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construction, each of these elements are numbers between 0 and 1, bounds included, 
as required for mole fractions. Also, the sum over the elements of any colunm of the 

x^^^ matrix is equal to 1. By construction, the matrix also relates the mixture 

molar densities of the phases and the components as described by equation (S), that is, 

5 m^^^^f^^rf^^ (35) 

The method ensures that mass is balanced because of the properties of the constructed 

and x^^^ matrices whose elements exhibit the properties described in equations 

(6) and (3), respectively. 

At block 1 12, the derivatives of the phase properties with respect to the 

10 mixture molar densities of the flowing components ^ needed for the iterative 

resolution of the system of equations that includes equations of the type of equation 
(23), are easily evaluated using the chain rule, which is known to those skilled in the 
art. For example, considering the intrinsic phase molar densities , for any phase 

/ and any flowing component 



15 



w J3» arnSa M> '■'^ 



where the derivatives with respect to the mixture molar densities of the EOS 
components, that is., d^jj dm^^^ , are provided by the phase equilibrium 

calculations. The fluid flow equations can then be solved. 

From the foregoing teaching of this invention, one skilled in the art may be 
20 inclined to peform flow calculations in the EOS domain with the mverse 

transformation replaced by the use of the gen^alized inverse, T^^'^ , of the 

transformation matrix 7^^^ , applied to the vector of the EOS conq)onents mixture 

molar densities "m^'^ obtained after the compositional miscible fluid flow has been 
performed: 

25 "wW=r(^).-wW=[r(^)'.^^)^].r(^f. -^w (s?) 



wo 02/057901 



21 



PCT/US02/01041 



However, because an EOS component can be present in several flowing components, 
this operation may not ensure that all the elements of the new vector of the mixture 
molar densities of the flowing components would be non-negative as required by the 
physics. This situation will be illustrated in the following simplified example. 
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Example calculation 

The following simplified example is provided to illustrate one aspect of the 
present invention. Only four flowing conqponents are considered: dry gas, wet gas, oil 
and water (i.e., JV^^^= 4). The normalized vector of the mixture molar densities of the 

5 flowmg components m^^^is n[^^ = (0. 10, 0.50,0.20,0.20)^ , respectively for the dry 
gas, the wet gas, the oil and the water The dry gas is simply methane (C,) . The wet 
gas is assumed to be composed of methane, propane (C3) and n-butane (C4) . The oil 
is assumed to be simply a combination of the paraflBns /jQo and wQg . The water 
component is common to both domains. The number of EOS components is thra 

li) N^/^ = 6 so that the transformation matrix 7^^^ is 6x4, chosen to be: 





0.8510 


0 


0> 


0 


0.1190 


0 


0 


0 


0.0300 


0 


0 


0 


0 


0.4440 


0 


0 


0 


0.5560 


0 


.0 


0 


0 


1^ 



It is easy to verify that the elements of the above matrix verify the required 
properties listed in equations (25) and (27). Application of equation (24) gives 
the normalized vector of the mixture EOS component molar densities: 

C, C3 nC^ nC,o nC,, H,0 

^(•)^j(*'),^(/) 52550 0.05950 0.0150 0.00880 0.1120 0.20/ 
As required by equation , mass balance is preserved as: 

Phase equilibrium is performed using the Peng-Robinson equation of state for 
the petroleum components at a pressure of 3000 psia (20,685 kPa) and a temperature 
20 of 212 T? (100 X). The sohibility of the petroleum components in the aqueous phase 
is governed by Henry's law. The calculations are performed in double precision and 
are ^ven with the scientific notations ZB+xx or E-xx that implies that the number Z 
has to be multiplied by 10 to the power +xx or -xx. For reason of presmtation, not all 
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the digits will be given thereafter. However, errors, if any, with respea to the 
expeded properties are given. 

The equilibrium corresponds to the presence of three phases: aqueous, liquid 

and vapor with a normalized mixture phase molar density tf'^ equal to 

5 Aqueous Liquid Vapor 

(19.411312951B.02 51,467908926E-02 29. 120778 123E-02)' 

and a non-square 6x3 composition matrix which is gjiven in Table 1. The sum 
over the elements for each of the the three columns is equal to 1 with no error. 

The non-square 3x6 matrix j^^'^ whose elements^ calculated using eqiiation 

10 (10), represent the distribution of the EOS conq)onents among the phases are given in 
Table 2. 

For a given component, the sum of the elements over all three phases is equal to 1 
with an error that never exceed 5 x 10"" % (see Table 3). 

The elements of the non-square matrix 3x4 y^^^ , calculated usmg equation 

15 ( 32 ) , are given in Table 4. All the elements are non-negative and less than 1 . As 
shown in Table S, for any flowing component, the error with respect to 1 of the sum 
of the elements over the phases is less than 44x10^" % , i.e., it is zero. Comparing 
Table 5 and Table 3, one may notice that the errors are the same for the dry gas and 
the water in both domains, as expected. 

20 The elements of the non-square 4x3 matrix x^^^ , computed using equation ( 9 ) 

, are listed in Table 6. All the numbers are non-negative and less than 1. The relative 
error compared to 1 of the sum of the elements over the components for each phase is 
-1 1. 102230246E-17% for the aqueous phase and zero for both the liquid and vapor 
phases. 

25 Another check can be made by considering the matrix product C^^^ =y^^^ • j^^^ 

(equation { 12 ) ) and its inverse, C^^^ . The elements of the 3x3 matrix C are given 
in Table 8. 

The sum over the elements of each cohunn is equal to 1 with a relative error of 
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(44.098058538E-14. 13.9888101 lOE-14, 46.629367034E-15)% from the left to the 
right, i.e., it is negligible. Using the vector of the mixture phase velocity. Hie product 
(<f) . ^(«) yields, as ejq)ected, the vector rf'^ with a relative error for the aqueous 
phase of 44.054170209E-14% for the aqueous phase, -13.61 1408413E-14% for the 
5 liquid phase and -53. 184056854B-14% for the v^r phase. These errors are 
negli^ble. 

Table 8 gives the dements of the 3x3 inverse matrix . As e}q)ected, the 
sum over the elements of each column is equal to 1, with again negligible relative 
errors of -45.996539910E-14»/o for the first column, 28.19966492E-14% for 

10 the second column and -19.140244945E-14% for the third one. 

For illustration of the usefulness of the present invention, Table 9 gives the 
elements of the generalized inverse of the composition matrix as indicated by 
equation (1 1). The mole fiiaction of the dry gas in the liquid phase is negative, and, as 
sudi, is physically unacceptable. Despite this problem, it is interesting to note that the 
15 sum over the phases for each of the flowing components is equal to 1 withui the 
precision of the calculations. 

Another way that could be argued to be considered is based on the 
transformation equation (24). Using the definition of the gena-alized inverse (see 
equation (1 1) in the case of the x matrix), one obtains: 

20 nl'^ = i^^-J'\ 

where T^^^ is the generalized inverse of the non-square matrix . Its elements 
aiB ven in Table 9. Usmg the rdation (5) in both domains along with the feet that 
the phases are identical in both domains, identification of the coefiBcients of the 
el^nents of the mbdure phase molar density rj yields: 

The elements of the non-square 4x3 matrix are ^ven in Table 10. 
Compared to Table 6, the mole fraction of the dry gas in the liquid phase is negative. 
This example shows that this procedure is only approximate. The problem for the dry 
30 gas is related to the feet that it is present in two of the flowing components. 
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WMe much of fhe foregoing description is directed toward oil field basin and reservoir 
simulation, it is understood that the disclosed method will also spply to models and gimnlatt/yng in a 
wide variety of sciences. Apeisonsldn6dinfiieart»particiilar^onehann^ 
teachings of this patent, wiU recognizB many modification 
5 disclosed above. As discussed above, fhe qtecificaUy disclosed embodinie^ 

be used to limit or restrict the scope of the invention, which is to be determinedby the claims bdow 
and tiieir equivaloits. 



wo 02/057901 



30 



PCT/US02/01041 



APPENDIX 

Superscripts 

: indicates that the variable belongs to the EOS (Equation of State) domain. 
5 : indicates that the variable belongs to the flowing components domain. 
When no superscript is present, the relation is valid in any domain. 
* : indicates that the matrix is a generalized inverse. 
^ : indicates the transpose of a matrix. 
: indicates the inverse of a matrix. 

10 

Subscripts 

A lower case Greek letter for example) indicates a component. 
An upper case Latin letter (j for example) indicates a fluid phase. 

15 Mathematical symbols 

V means "for all". 
€ means 'belongs to'. 

[1, J^T] means 'interval 1 to J^, bounds included'. 
A vector is indicated by an underlined letter (m for example). 
20 A matrix is indicated by a double underlined letter (x for example). 

^ X means sum of ejqpression from a=l to . 

The Einstein convention is not used C e., no sum over repeated indices), 
div indicates the divergence operator (with respect to Eulerian coordinates), 
gmd indicates the gradient op^ator (with respect to Eulerian coordinates). 

25 A matrix product is indicated by the symbol ' \ 
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Greek letters 

Sjj : indicates the Kronecker symbol (=0 for /^t J , =1 for /= J ). 
fj : Np vector ofthe mixture molar densities of the phases. 
rij : mixture molar density of fluid phase J. 
5 jiij : viscosity of fluid phase J. 

4> : intrinsic molar density of fluid phase J. 
Pj lintrinsicmassdensity of fluid phase J. 
^ : porosity of material. 

10 Latin letters 

c : N^xNg matrix equal to the product x- y* . 

c : N^xN^ matrix equal to the product X- j^. 

C* :NpxNp matrixequaltotheproduct /-x. 

C : Np xNp matrix equal to the product y- x. 

IS d^j : mixture molar density of component a in fluid'phase J. 

g : acceleration of gravity. 
: Np^Np identity matrix. 

UTj : relative permeability offluid phase 7. 

K : absolute permeability of a material 
20 m : total mixture molar density of a system. 

m : vector ofthe mixture molar density of the components. 

: mixture molar density of conq)onent a . 
: number of components. 
Np : number of phases. 
25 P : reference pressure. 

P^j : capillary pressure of fluid phase J. 
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^ : source mixture molar density rate for component 
Sj : saturatioii of fluid phase J. 

vector whose elements are the particle velocities of the components 

5 : vector of the particle velocity of fluid phase 7. 

F^^^ : vector of the particle velocity of the solid phase. 
t : time. 

: JV^'^ matrix of the mapping from the floi^g conq)onents domain to the 

EOS components domain. 
10 x^j ; mole fraction of component a in fluid phase J, (composition of the phase). 

X : xNp matrix whose elements are the x^ . Each cohmm gives the composition 

of 

a phase. 

yj^ : fluid phase (J) mole fraction of component a , 
IS y : NpxN^ matrix whose elements are the yj„ . Each column gives the distribution 

of 

a component among the phases. 
z : depth (Eulerian coordinate). 
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What is daimed is: 

1. A method of estimatmg one or more properties of a multi-component fluid 
contained in a physical system containing two or more phases, comprising the 
steps of 

S (a) equating the physical system in at least one dimension to a multiplicity 

of cells; 

(b) for each cell, characterizing the multi-component fluid in the cell using a 
property of a first set of components and representing the 
characterization by a first vector; 
10 (c) converting the first vector to a second vector usmg a transformation 

matrix, the transformatin matrix being indicative of the distribution of 
the first set of components and the second vector being representative of 
a property of a second set of components greater in number than the first 
set of components; 

15 (d) using the second vector to determine the number and properties of 

phases present in each cell; 

(e) determining the elements of a first matrix that expresses distribution of 
the second vectors among the phases; 

(f) determining the elements of a second matrix that expresses how the first 
20 vectors are distributed among the phases using the transformation matrix 

and the first matrix; 

(g) determining a third matrix that expresses the composition of the phases 
in terms of the first vectors; and 

(h) using the third matrix to perform fluid flow calculations to estimate one 
25 or more properties of the multi-conq)onent fluid. 

2. A method for emulating fluid movement and phase change behavior, wherein 
the method comprises: 

(a) determining for each of a plurality of discrete cells a vector of first 
30 domain component concentrations; 

(b) transforming the first domain component vectors into vectors of second 
domain component concentrations; 
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10 



(c) applying to fhe second domain component vectors, equations of state to 
detennine second domain component distribution matrices indicative of 
a distribution of the second domain components among eacht of the 
phases; 

(d) converting the second domain component distribution matrices to first 
domain component distribution matrices; and 

(e) ^plying the first domain component distribution matrices to fluid flow 
equations to determine for each of the plurality of discrete cells an 
updated vector of first domain con^onent concentrations. 

The method of claim 2, wherein the converting indudes multiplying the 
second domain component distribution matrix with a transform matrix. 



4. The method of claim 2, further comprising: 

15 displaying the concentrations of a desired component in a given cell at a given 

time. 

5. A nonlinear process simulator which comprises: 
(a) a display monitor; 

20 (b) a processor coupled to the display monitor and configured to display 

time evolution of the process; and 
(c) a memory configured to store soflware for access by the processor, 
wh^ein the software includes: 

a second domain module which configures the processor to ddiermine a 
25 distribution of second domain components among a phirality of phases, 

wherein the second domain conq)onent distribution is expressible in 
terms of a second domain component distribution matrix having 
elements that specify what fraction of each second domam component is 
found in a ^ven phase; and 
30 a conversion module which configures the processor to convert the 

second domain component distribution matrix into a first domain 
component distribution matrix by multiplying the second domain 
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component distribution matrix with a transform matrix to obtain a first 
domain component distribution matrix, wherein the second domain 
components are expressible as a product of the first domain components 
and the transform matrix. 

5 

6. The simulator of claim 5 wherein the software fiirther includes: 

(a) a first domain module which configures the processor to detemune for 
each of a plurality of discrete cells a vector of first domain component 
concentrations; 

10 (b) a transform module which configures the processor to transform the 

first domain component concentration vectors into second domam 
component concentration vectors^ wherein the second domain module 
operates on the second domain component concentration vectors to 
detenfiine the second domain component distribution matrix 

15 

7. The simulator of claim 6 wherein the software finther includes a display 
module that configures the processor to display concentrations of a desired 
component in a given cell at a given time. 

20 8. The simulator of claim 6, wherein the first domain module performs 

simulation of fluid flow, and the second domaui module performs simulation 
of fluid phase behavior. 



9. The simulator of claim 8, wherein the software configures the processor to 
25 simulate development of petrolwm reservoirs by repetition of actions that 

include simulating fluid flow in the first domain, transforming first domain 
component concentrations into second domain component concentrations, 
simulating fluid phase behavior in the second domain, and converting the 
second domam component distribution matrix into the first domain component 
3b distribution matrix. 
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10. The simulator of claim 8, wherein the software configures the processor to 
simulate movement of petroleum in petroleum reservoirs by repetition of 
actions that include simulating fluid flow in the &st domain, transforming first 
domain component concentrations into second domain component 

5 concentrations, simulating phase behavior in the second domain, and 

converting the second domain conqponent distribution matrix into the first 
domain component distribution matrix. 

11. An information carrier medium configured to provide a processor with a 

10 program, wherdn when the processor executes the program the processor is 

configured to: 

(a) determine a distribution of second domain components among a plurality 
of phases, wh^ein the distribution is e?q)ressible as a second domain 
component distribution matrix having elements that specify what 

15 fraction of each second domain component has a given phase; and 

(b) convert the distribution of second domain components among the phases 
to a distribution of first domain components among the phases, wherein 
the second domain components are expressible as a product of the first 
domain components and a transform matrix, and wherein the conversion 

20 of the distribution includes multiplying the second domain component 

. distribution matrix with the transform matrix to obtain a first domain 
component distribution matrix. 

12. The carrier medium of claim 1 1, wherein the carrier medium is a computer-- 
25 readable information storage medium. 

13 . The carrier medium of claim 1 1, wherein the carria: medium is an information 
transmission medium. 



30 14. 



The carrier medium of claim 1 1, wherein the program fiirther configures the 
processor to: 
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(a) evaluate equations expressed in terms of a first domain to detmnine first 
domain component values; 

(b) transform the first domain component values into second domain 
component values; and 

5 (c) evaluate equations expressed interna of a second domain to determine 

the distribution of second domain components. 



15. The carrier medium of clmm 14, wherein the first domain equations include 
fluid flow continuity equations, and wherein the second domain equations 
10 include phase equilibrium and flash equations. 



16. The carrier medium of claim 15, wherein the program configures the processor 
to simulate development of petroleum reservoirs by repetition of actions that 
include evaluating fluid flow equations in the first domain, transforming first 
15 domain component values into second domain component values, evaluating 

phase equilibrium and flash equations in the second domain, and converting 
the second domain component distribution into the first domain component 
distribution. 



20 17. The carrier medium of claim 15, wherein the program configures the processor 
to simulate movement of petroleum in petroleum reservours by repetition of 
actions that include evaluating fluid flow eqxiations in the first domain, 
transforming first domain component values into second domain component 
values, evaluating phase equilibrium and flash equations in the second 

25 domain, and converting the second domain component distribution into the 

first domain component distribution. 
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